% Tables 1, 2 and 8 in "Commodity Prices, Convenience Yields, and
% Inflation," March 2011

clear all;

load cyj.dat;                   % convenience yields on 23 commodities
load qj.dat;                    % HP-detrended real commodity prices
load us_intrate_monthly.dat;    % US 3-month T-bill rate
load us_cpi_monthly.dat;        % US CPI indeces (all items, lfe etc.)
load imf_index_monthly.dat;     % IMF aggregate commodity price index
load usd_monthly_ave.dat;       % USD trade-weighted index
load us_unempl.dat;             % US unemployment rate

i=us_intrate_monthly;
P=us_cpi_monthly;
SA=imf_index_monthly;
er=usd_monthly_ave;
ur=us_unempl;
nn=rows(P);

rhat=2;         % number of principal components
nwl=-1;         % number of lags in Newey-West HAC estimator
bsize=4;        % bootstrap block size
B=4999;         % number of bootstrap replications

[ytu,gapu]=hptrend(ur,14400);                       % HP-detrending
[ehat,Fhat,lamhat,ve]=pc(standard(cyj),rhat);       % PC estimation
[ehat2,Fhat2,lamhat2,ve2]=pc(standard(qj),rhat);    % PC estimation

for inflmeas=1:6;     % inflation measure (all items, lfe etc.)
    
if inflmeas==1;
    hh=[1 3 6 12]';   % forecast horizon
else
    hh=[1 3]';
end

for hi=1:rows(hh);
h=hh(hi);

ds_h=100*(log(SA(3+h:nn,1))-log(SA(3:nn-h,1)));
ds_1=100*(log(SA(3:nn-h,1))-log(SA(2:nn-h-1,1)));
ds_2=100*(log(SA(2:nn-h-1,1))-log(SA(1:nn-h-2,1)));
dx=100*(log(er(3:nn-h,1))-log(er(2:nn-h-1,1)));
ir=i(3:nn-h,1);
gap=gapu(3:nn-h,1);

infl_h=100*(log(P(3+h:nn,inflmeas))-log(P(3:nn-h,inflmeas)));
infl_1=100*(log(P(3:nn-h,inflmeas))-log(P(2:nn-h-1,inflmeas)));
infl_2=100*(log(P(2:nn-h-1,inflmeas))-log(P(1:nn-h-2,inflmeas)));
p1_cy=Fhat(3:nn-h,1);
p2_cy=Fhat(3:nn-h,2);
p1_rcpd=Fhat2(3:nn-h,1);
p2_rcpd=Fhat2(3:nn-h,2);
roil1=qj(3:nn-h,22);
roil2=qj(2:nn-h-1,22);

% Commodity price equation
x1=[ones(rows(p1_cy),1) p1_cy p2_cy ir dx ds_1 ds_2];
res1=nwest(ds_h,x1,0);
b1=res1.beta;
r1=res1.resid;
R2_1=res1.rbar;
G1=nw((r1*ones(1,cols(x1))).*x1,nwl);
V1=inv(x1'*x1)*G1*inv(x1'*x1);
se1 = sqrt(diag(V1));
ts1 = b1./se1;

% Baseline inflation equation
x=[ones(rows(p1_cy),1) p1_cy p2_cy p1_rcpd p2_rcpd infl_1 infl_2];
res=nwest(infl_h,x,0);
b=res.beta;
r=res.resid;
R2=res.rbar;
G=nw((r*ones(1,cols(x))).*x,nwl);
V=inv(x'*x)*G*inv(x'*x);
se = sqrt(diag(V));
ts = b./se;

% Augmented inflation equation
if inflmeas==1;
    x2=[ones(rows(p1_cy),1) p1_cy p2_cy p1_rcpd p2_rcpd infl_1 infl_2 roil1 roil2 ir dx gap];
else;
    x2=[ones(rows(p1_cy),1) p1_cy p2_cy p1_rcpd p2_rcpd roil1 roil2 infl_1 infl_2];
end;
res2=nwest(infl_h,x2,0);
b2=res2.beta;
r2=res2.resid;
R2_2=res2.rbar;
G2=nw((r2*ones(1,cols(x2))).*x2,nwl);
V2=inv(x2'*x2)*G2*inv(x2'*x2);
se2 = sqrt(diag(V2));
ts2 = b2./se2;

% Bootstrap procedure
bmatrix=[cyj(3:nn-h,:) qj(3:nn-h,:) roil1 roil2 infl_h infl_1 infl_2 ds_h ir dx gap ds_1 ds_2 p1_cy p2_cy p1_rcpd p2_rcpd];
gg=floor(rows(bmatrix)./bsize);
bmatrix=bmatrix(1:gg*bsize,:);
betab=[]; betab1=[]; betab2=[];
tstatb=[]; tstatb1=[]; tstatb2=[];
for j=1:B;
    indexx=floor(1+gg*unifrnd(0,1,1,gg))';
    ind=blockind(indexx,bsize);
    ind=ind(1:rows(bmatrix));
    matrixb=bmatrix(ind,:);
    cyjb=matrixb(:,1:23);
    qjb=matrixb(:,24:46);
    roilb1=matrixb(:,47);
    roilb2=matrixb(:,48);
    inflb_h=matrixb(:,49);
    inflb_1=matrixb(:,50);
    inflb_2=matrixb(:,51);
    dsb_h=matrixb(:,52);
    irb=matrixb(:,53);
    dxb=matrixb(:,54);
    gapb=matrixb(:,55);
    dsb_1=matrixb(:,56);
    dsb_2=matrixb(:,57);
    p1_cyt=matrixb(:,58);
    p2_cyt=matrixb(:,59);
    p1_rcpdt=matrixb(:,60);
    p2_rcpdt=matrixb(:,61);
    
    [ehatb,Fhatb,lamhatb,veb]=pc(standard(cyjb),rhat);
    [ehatb2,Fhatb2,lamhatb2,veb2]=pc(standard(qjb),rhat);
    c1=corr([Fhatb(:,1) p1_cyt]);
    c2=corr([Fhatb(:,2) p2_cyt]);
    c3=corr([Fhatb2(:,1) p1_rcpdt]);
    c4=corr([Fhatb2(:,2) p2_rcpdt]);
    p1_cyb=sign(c1(1,2))*Fhatb(:,1);
    p2_cyb=sign(c2(1,2))*Fhatb(:,2);
    p1_rcpdb=sign(c3(1,2))*Fhatb2(:,1);
    p2_rcpdb=sign(c4(1,2))*Fhatb2(:,2);
    
    x1b=[ones(rows(p1_cyb),1) p1_cyb p2_cyb irb dxb dsb_1 dsb_2];
    [b1b,bint1b,r1b,rint1b,stats1b] = regress(dsb_h,x1b);
    G1b=nw((r1b*ones(1,cols(x1b))).*x1b,nwl);
    V1b=inv(x1b'*x1b)*G1b*inv(x1b'*x1b);
    se1b = sqrt(diag(V1b)); 
    betab1(j,:)=(b1b-b1)';
    tstatb1(j,:)=((b1b-b1)./se1b)';
    
    xb=[ones(rows(p1_cyb),1) p1_cyb p2_cyb p1_rcpdb p2_rcpdb inflb_1 inflb_2];
    [bb,bintb,rb,rintb,statsb] = regress(inflb_h,xb);
    Gb=nw((rb*ones(1,cols(xb))).*xb,nwl);
    Vb=inv(xb'*xb)*Gb*inv(xb'*xb);
    seb = sqrt(diag(Vb));
    betab(j,:)=(bb-b)';
    tstatb(j,:)=((bb-b)./seb)';
    
    if inflmeas==1;
        x2b=[ones(rows(p1_cyb),1) p1_cyb p2_cyb p1_rcpdb p2_rcpdb inflb_1 inflb_2 roilb1 roilb2 irb dxb gapb];
    else;
        x2b=[ones(rows(p1_cyb),1) p1_cyb p2_cyb p1_rcpdb p2_rcpdb roilb1 roilb2 inflb_1 inflb_2];
    end;
    [b2b,bint2b,r2b,rint2b,stats2b] = regress(inflb_h,x2b);
    G2b=nw((r2b*ones(1,cols(x2b))).*x2b,nwl);
    V2b=inv(x2b'*x2b)*G2b*inv(x2b'*x2b);
    se2b = sqrt(diag(V2b));
    betab2(j,:)=(b2b-b2)';
    tstatb2(j,:)=((b2b-b2)./se2b)';

end;

fprintf(' inflmeas = %6.4f\n',inflmeas);
fprintf(' h = %6.4f\n',h);
if inflmeas==1;
fprintf('\n estimation results for commodity price model\n')
fprintf('    coeff.     s.e.    boot s.e.  90%% bootstrap CI\n')
% equal-tailed percentile-t bootstrap
disp([b1 se1 std(betab1)' b1-se1.*quantile(tstatb1,0.95)' b1-se1.*quantile(tstatb1,0.05)'])
% Hall's percentile bootstrap
%disp([b1 se1 std(betab1)' b1-quantile(betab1,0.95)' b1-quantile(betab1,0.05)'])
fprintf(' Rbar^2 = %6.4f\n',R2_1);
fprintf('\n estimation results for baseline inflation model\n')
fprintf('    coeff.     s.e.    boot s.e.  90%% bootstrap CI\n')
% equal-tailed percentile-t bootstrap
disp([b se std(betab)' b-se.*quantile(tstatb,0.95)' b-se.*quantile(tstatb,0.05)'])
% Hall's percentile bootstrap
%disp([b se std(betab)' b-quantile(betab,0.95)' b-quantile(betab,0.05)'])
fprintf(' Rbar^2 = %6.4f\n',R2);
end
fprintf('\n estimation results for augmented inflation model\n')
fprintf('    coeff.     s.e.    boot s.e.  90%% bootstrap CI\n')
% equal-tailed percentile-t bootstrap
disp([b2 se2 std(betab2)' b2-se2.*quantile(tstatb2,0.95)' b2-se2.*quantile(tstatb2,0.05)'])
% Hall's percentile bootstrap
%disp([b2 se2 std(betab2)' b2-quantile(betab2,0.95)' b2-quantile(betab2,0.05)'])
fprintf(' Rbar^2 = %6.4f\n',R2_2);

end

end